function h = jcb_plot_index(filename,draw_bowtie,draw_crystal,colorname,savefigure,fileformat)

%% initial setup
crystal_middle_size=35;
num_theta_points=100;

if nargin<6
    fileformat='pdf';
end
if nargin<5
    savefigure=false;
end
if nargin<4
    colorname='red';
end
if nargin<3
    draw_crystal=true;
end
if nargin<2
    draw_bowtie=true;
end

if nargin<1
    % get user to pick files
    [filelist,pathname] = jcb_pick_files('b*index_data.txt', ...
    'Please select files to create peak list for',true,pwd);
else
   if strcmpi(filename,'all')
        filename='b';
    end
    [filelist,pathname] = jcb_select_files(filename,'index');
end
filelist=filelist';

%%
for ii=1:size(filelist,1)
    %% load file & extract data
    filename=strtrim(filelist{ii});
    disp(['creating figure for ' filename]);
    profile = load(fullfile(pathname,filename));
    lambda=profile(1,1);
    x=profile(2:end,1);
    y=profile(1,2:end);
    offsetx=diff(x)/2;
    offsety=diff(y)/2;
    offsetx=offsetx([1:end end]);
    offsety=offsety([1:end end]);
    profile=profile(2:end,2:end);
    
    %% plot profile
    h = figure(63256);
    clf;
    jcb_contourplot(x-offsetx,y-offsety,profile);
    caxis([0, max(caxis)]);
    colormap(jcb_colormap([],4));
    title(['refractive index profile for ' strrep(filename,'_','\_')]);
    shading faceted


    %% draw crystal:
    if draw_crystal
        dose=regexp(filename,'dose(?<btl>\d+)_(?<btw>\d+)_(?<gap>\d+)_','names');
        if isempty(dose)
            dose=regexp(filename,'dose(?<btl>\d+)','names');
            dose.btl=dtmp.btl;
            dose.btw=dtmp.btl;
            dose.gap=num2str(20+(80-str2double(dtmp.btl)));
        end
        fileinfo.btl=str2double(dose.btl);
        fileinfo.btw=str2double(dose.btw);
        fileinfo.gap=str2double(dose.gap);
        mesh=regexp(filename,'m(?<mesh>\d+)','names');
        fileinfo.mesh=str2double(mesh.mesh);
        loc=regexp(filename,'loc(?<loc>\d+)sweep_d(?<d>\d+)','names');
        if isempty(loc)
            loc=regexp(filename,'loc(?<loc>\d+)','names');
        end
        fileinfo.crystal_x=0;
        fileinfo.crystal_angle=45/180*pi;
        switch str2double(loc.loc)
            case 1
                if length(fields(loc))>1
                    fileinfo.crystal_y=str2double(loc.d);
                else
                    fileinfo.crystal_y=10;
                end
            case 2
                fileinfo.crystal_y=20;
            case 3
                fileinfo.crystal_y=40;
            case 4
                fileinfo.crystal_y=108;
            case 5
                fileinfo.crystal_x=75;
                if length(fields(loc))>1
                    fileinfo.crystal_y=str2double(loc.d);
                    fileinfo.crystal_angle=atan(0.5);
                else
                    fileinfo.crystal_y=71;
                end
            case 6
                fileinfo.crystal_y=[];  % as no crystal used
        end

        
         % check profile height
        fid=fopen(fullfile(pathname,filename),'r');
        fileinfo.height=fscanf(fid,'%c',300);
        fileinfo.height=regexp(fileinfo.height,'height z=(\d+.?\d+)','tokens');
        fileinfo.height=str2double(fileinfo.height{1});
        fclose(fid);
        
        % calculate the crystal cross-section size at the given height
        cryst_size=2*(10+(35-fileinfo.height)*(crystal_middle_size-20)/2/17.5);

        % define crystal coordinates, centred at origin
        coords = 0.5*[-cryst_size, -cryst_size;
                  -cryst_size,  cryst_size;
                   cryst_size,  cryst_size;
                   cryst_size, -cryst_size;
                  -cryst_size, -cryst_size]; 

        % rotate crystal
        for kk=1:size(coords,1)
            coords(kk,:) = ([cos(fileinfo.crystal_angle), -sin(fileinfo.crystal_angle);
                            sin(fileinfo.crystal_angle),  cos(fileinfo.crystal_angle)]*coords(kk,:)')';
        end
        % shift coords to the correct position for the crystal
        coords = coords + ones(size(coords,1),1)*[fileinfo.crystal_x fileinfo.crystal_y];

        if ~isempty(fileinfo.crystal_y)
            line( coords(:,1), coords(:,2),'color',colorname);
        end
    end

    %% draw bowtie
    if draw_bowtie
              
        p  = 2*atan(fileinfo.btw/2/fileinfo.btl);
        q  = (pi-p)/2;
        Rc = 9;
        L  = fileinfo.btl;
        W  = fileinfo.btw/2;
        L=L+Rc*(1/sin(p/2)-1);
        W=W+Rc*(1/tan(q/2)-1);


        inset=Rc/sin(p/2);
        Wa = Rc*cos(p/2);
        Ha = Rc*sin(p/2);

        e  = Rc/tan(q/2);
        f  = e*cos(q);
        g  = e*sin(q);

        lineX = [ [inset-Ha, L-g];
                  [L, L]
                  [L-g, inset-Ha]] - (inset-Rc)+fileinfo.gap/2;
        lineY = [[Wa, W-f];
                 [W-e, -W+e]
                -[W-f, Wa]];

        circleX = zeros(3,num_theta_points);
        circleY = zeros(3,num_theta_points);
        theta=linspace(2*pi-(pi/2+p/2),pi/2+p/2,num_theta_points);
        circleX(1,:) = Rc*cos(theta)+Rc+fileinfo.gap/2;
        circleY(1,:) = Rc*sin(theta);
        theta=linspace(pi/2+p/2,0,num_theta_points);
        circleX(2,:) = Rc*cos(theta)+L-inset+fileinfo.gap/2;
        circleY(2,:) = Rc*sin(theta)+W-e;
        theta=linspace(0,-pi/2-p/2,num_theta_points);
        circleX(3,:) = Rc*cos(theta)+L-inset+fileinfo.gap/2;
        circleY(3,:) = Rc*sin(theta)-W+e;
        fullX=[circleX(1,:) lineX(1,:),circleX(2,:) lineX(2,:),circleX(3,:) lineX(3,:)];
        fullY=[circleY(1,:) lineY(1,:),circleY(2,:) lineY(2,:),circleY(3,:) lineY(3,:)];           
        line(fullX,fullY,'color',colorname);
        line(-fullX,fullY,'color',colorname);

    end
    
    %% save figure to file
    if savefigure
        print(['-d' fileformat],regexprep(filename,'_data\.txt',''));
    end
end

if nargout<1
    clear h
end
end